The invention relates to a method of predicting acoustic wave performance in sediments. In particular, the invention relates to a compacted shale model and applications of the model to generate theoretical sonic logs useful in seismic studies. The model predicts acoustic velocity with depth (sonic log) and allows procedures such as sonic log editing and quality control.
Acoustic wave propagation in sedimentary rock sequences, the subject of sonic logs, is of fundamental interest in petroleum exploration. It provides a key linkage between geophysical data acquisition and interpretation, and the rock properties which are of basic interest to geologists. Compression (and more recently shear) wave data are commonly acquired during borehole logging operations. These data are subsequently used in and are often critical to interpretation of rock properties, reservoir analyses, seismic interpretation and basin analyses.
The linkage between rock properties measured in boreholes and the interpretation of similar properties from seismic data is provided by measurements recorded in-situ by logging tools, and by detailed laboratory measurements on rock samples recovered during drilling. The properties of greatest interest in this process are acoustic performance (of both compression and shear waves) and bulk density. The borehole environment and logging process often adversely affect acquisition of good quality acoustic log data over substantial intervals of section, resulting in poor ties of well data to seismic, and inferior quality acoustic velocity data.
The quality of the borehole log data is often affected by petrophysical properties (fractures, compaction, hydrocarbon content), borehole environmental factors (mud properties, borehole surface conditions) and acquisition parameters (logging speed, signal generation and detection techniques). The data so acquired are calibrated by comparing the integrated borehole signal with independently measured interval velocity data (check-shot data). Misfit between integrated log data and check-shot data primarily arises because noise is commonly incorporated in the log signal, and because of different acoustic frequencies employed in the two techniques (Ward, R. W. and Hewitt, M. R., 1977, Monofrequency borehole travel time survey: Geophysics, 42, 1137-1145). Goetz et al (1979, An investigation into discrepancies between sonic log and seismic check-shot velocities: Australian Petroleum Exploration Association. J.. 19, 2, 131-141) provided a complete discussion of error sources of the two processes.
Corrections are applied to the borehole data to force a fit to the check-shot data, and the emergent calibrated sonic log is then used as input to further studies, particularly seismic modelling. The tie of well data to seismic, and the interpretation of seismic character of sedimentary packages (seismic stratigraphy) is fundamental to interpretation of basic structural evolution, the history of deposition and present geometry.
Sonic signal degradation, particularly in near-surface and less compacted rocks, often leads to substantial editing being required before the integrated signal agrees acceptably with the check-shot data. Lack of a suitable technique (both in terms of physical modelling and operational efficiency) for systematic noise removal and editorial replacement of intervals of suspect data has hitherto resulted in linear interpolation being the most commonly used method of noise removal from the sonic log.
Under normal circumstances, in generally subsiding depositional basins, progressively increasing overburden load due to increasing depth of burial, results in sequence compaction, with porosity reduction, increased bulk density and improved acoustic propagation efficiency. Observation of sonic logs clearly shows a general increase in acoustic p-wave velocities of propagation with increasing depth of burial (Telford, W. M., Geldart, L. P., and Sheriff, R. E., 1990, Applied Geophysics, (Second Edition), Cambridge University Press). Exceptions exist, and are primarily lithology-dependent. Several recent papers have used various compaction models to quantify these changes, and to use the data for studies of basin evolution.
An exponential decay model for the density-depth function was proposed by Stegena, L. (1964, The structure of the earth""s crust in Hungary, Acta Geologica, Budapest 8, 413-431), and Korvin, Gabor (1984, Shale compaction and statistical physics: Geophysical Journal of the Royal Astronomical Society, v. 78, p. 35-50) developed a mathematical proof of the exponential decay model for shale compaction. This model has not been widely used (Japsen, Peter, 1998, Regional velocity-depth anomalies, North Sea Chalk: a record of overpressure and Neogene uplift and erosion: AAPG Bulletin, v82, No 11, p. 2031-2074, Heasler, Henry P., and Kharitonova, Natalya A., 1996, Analysis of sonic well logs applied to erosion estimates in the Bighorn Basic, Wyoming: AAPG Bulletin, v. 80, No. 5. p. 630-646). Difficulties in the use of such a model arise from the often complex mix of lithologies and absence of observational data for most lithologies other than shale.
Japsen uses a segmented linear model for North Sea Chalk. Gassman, F., (1951, Ueber die elastizitat poroser medien: Natur. Ges. Zurich, Vierteljahrssch. V. 96, p. 1-23) introduced a physical model for compressional wave velocity in porous rocks, and this has been recently applied to quantify variations in sonic p-wave performance in sandstone reservoirs (Alberty, Mark, 1996, The influence of the borehole environment upon compressional sonic logs: The Log Analyst, v. 37, p. 30-44).
It is an object of the invention to use data acquired in association with boreholes in an improved manner by means of mathematics based processing to generate synthetic sonic logs. It is a particular object of the invention to provide methods by which to overcome problems such as the defects in actual logs, which logs are often compromised by borehole engineering, environmental difficulties, and by operational considerations. More particularly, the invention may provide an improvement over the considerable post-acquisition editing and re-calibration of the sonic signal which has been required so as to yield acceptable agreement between integrated sonic and check-shot measured interval travel times.
A study of the existing methods of editorial calibration, and a consideration of the underlying physical and mathematical processes, has led to a re-examination of the physics and mathematics of compactive modelling, and to the development of methods of use as outlined hereinbelow.
Because sedimentary rock burial history determines density, applying the invention allows interpretation of rock velocity in terms of burial history (the depth z in equation 4 below).
Systematic response of p-wave propagation efficiency to progressively increasing compaction is approached by initial consideration of the response of a relatively pure lithology. A marine shale with low total organic carbon content, buried progressively but sufficiently slowly that a normal pore pressure gradient is maintained, is first considered.
Compaction algorithm.
An exponential decay model is used below (after Korvin, 1984), wherein the shale density xcfx81 progressively changes with burial depth:
xcfx81(z)=xcfx81∞+(xcfx810xe2x88x92xcfx81∞)exe2x88x92kz xe2x80x83xe2x80x83(1) 
where xcfx81(z) is the bulk density at depth z, k is a compaction constant, ∞ and xcfx810 are respectively the bulk density at infinite depth and at the mudline and e is the exponential constant.
Boundary considerations yield clearly defined limits to the above. At the mud-line, as clastic debris (fragments of pre-existing rock) first accumulates, the newly deposited material will have a bulk density (xcfx810) similar to that of the water of deposition, that is, about 1.022 for seawater, and 1.000 for fresh water. Upon burial, initial consolidation is rapid, and bulk density increases accordingly. With increasing depth of burial, the density asymptotically approaches a limit (xcfx81∞) which is the upper limit of shale densityxe2x80x94about 2.7 gm/cc. ie.
xcfx81(z)=2.7+(xcfx810xe2x88x922.7)exe2x88x92kZ 
Solution for k is straightforward. This is best achieved by a least-squares best fit to data observed at a number of different depths. xcfx810 is established by consideration of depositional conditions, on geologic grounds. The constant k is the compaction constant, the larger the value of k, the more rapid the compaction (or vice versa). It is to some extent dependent on the geologic environment (particularly time and temperature). If depths are chosen in kilometers, and densities are expressed in normal units, according to Korvin, the value of k ranges from 0.28 kmxe2x88x921 to 1.46 kmxe2x88x921.
Similar compactive functional descriptors can be expected for other lithologies, though rock fabric variations and chemical stability considerations somewhat complicate the modelling process.
ACOUSTIC ALGORITHM.
An extension of the Korvin model to acoustic p-wave velocities, and the development of a practical technique to use this model to correct observed sonic data, follows.
Acoustic p-wave velocity profiles, measured by either move-out studies or by borehole acquisition techniques (acoustic logs and borehole check-shot surveys) show a similar trend of increasing velocity with depth. Because the rock sequences so measured are complex, velocity functions of various forms (though rarely exponential) have usually been developed to fit observed data, and provide the keys to time-depth conversion.
The velocity function for a normally compacted pure shale can be established as follows.
From an examination of wave theory (see, for example, Gorbachev, pp. 96-99)                               Vp          z                =                                                            [                                                      λ                    +                                          2                      ⁢                                              xe2x80x83                                            ⁢                      μ                                                                            ρ                    z                                                  ]                                            1                2                                      ⁢                          xe2x80x83                        ⁢            and            ⁢                          xe2x80x83                        ⁢                          Vs              z                                =                                    (                              μ                                  ρ                  z                                            )                                      1              2                                                          (        2        )            
where V is the velocity of sound and Vp is the velocity of a p wave, Vs is the velocity of shear wave, Vpz is the velocity of a p wave at depth z, Vsz is the velocity of an s wave at depth z, xcex is the Lamxc3xa9 constant related to bulk modulus and xcexc the shear modulus.
Since initial interest is in p-wave behaviour, substitute Vp(z) into the relationship of Korvin, equation 1 above:                     λ        z            +              2        ⁢                  μ          z                                    (                  Vp          ⁢                      (            z            )                          )            2        =                    (                              λ            ∞                    +                      2            ⁢                          μ              ∞                                      )                    Vp        ∞        2              +                  (                              λ            0                    +                      2            ⁢                          μ              0                                      )                    Vp        0        2              -                            (                                    λ              ∞                        +                          2              ⁢                              μ                ∞                                              )                          Vp          ∞          2                    ⁢              ⅇ                  -          kz                    
Now, for a single lithotype, xcex and xcexc at limits are constants, so the above reduces to:                                                         λ              z                        +                          2              ⁢                              μ                z                                                                        Vp              ⁡                              (                z                )                                      2                          =                              A                          Vp              ∞              2                                +                                    (                                                B                                      Vp                    0                    2                                                  -                                  A                                      Vp                    ∞                    2                                                              )                        ⁢                          ⅇ                              -                kz                                                                        (        3        )            
where A and B are constants.
Now for sonic velocities the interval transit time is   ITT  =            10      6        Vp  
(ITT is the standard petrophysical abbreviation for Interval Transit Time, the inverse of sonic velocity), so substituting in the above:
xcex94T(z)2*(xcexz+2xcexcz)=(A*xcex94T∞2)+(B*xcex94T02xe2x88x92A*xcex94T∞2)exe2x88x92kz xe2x80x83xe2x80x83(4) 
where xcex94T(z) is the p wave interval transit time (ITT) at depth Z, xcex94T0 is the ITT at the depositional surface, xcex94T∞ is the p-wave ITT at infinite depth.
Note that the constant k is precisely the same compaction constant as that derived for the density function above. Thus values of k derived from one set of observational data (eg. core data) can be applied in the other algorithm. Once k is determined, it can be used interchangeably in density and sonic algorithms. If the Lamxc3xa9 constants xcex and xcexc can be determined, equation (4) can be used to solve for xcex94T.
Within the single, progressively compacting lithotype model used above, sonic velocity is also observed to increase systematically with increasing depth of burial. This implies that the Lamxc3xa9 constants also change predictably and systematically. From boundary conditions solve       Vp    z    =            [                        λ          +                      2            ⁢                          xe2x80x83                        ⁢            μ                                    ρ          z                    ]              1      2      
for (xcex+2xcexc) at surface and at infinite (great) depth. At the mudline, xcex+2xcexc=2.32*106 and at great depth, xcex+2xcexc=70*106. The mixed term (xcex+2xcexc) is referred to as the elasticity function.
APPLICATION TO WELL DATA.
In a well section, if we have prior or independent knowledge of k (e.g. from core density observation), and an acceptable value for xcex94T in a clean shale lithology at a known depth, we can solve (4) above for (xcex+2xcexc) at that depth. We can gather such data for a number of depths, and develop by interpolation a continuous function for (xcex+2xcexc) with depth. This yields a method for computing xcex94T with depth.
Under normal circumstances, for clean shale, compaction is irreversible. Departures from the clean shale compaction trend arise from the inclusion of non-shale materials. The most common of these, resulting in deviations from the normal compactive trend, including apparent under-compaction, are organic carbon and water, or a mix of other lithologies.
Abrupt departure from the normal shale compaction trend towards higher density and velocity is generally due to uplift and erosion of part of the overburden. Inclusion of abnormally dense material (eg complex iron compounds which occur in oolites in the Evergreen Formation in the Surat Basin) can have a similar effect, but these events are generally short-term, and readily recognised. Where a stepwise velocity increase persists, an erosional break can be inferred, and may be interpreted in quantitative terms.
The above argument is applied to a single, progressively compacting shale lithology. The model presumes that the changes in shale density are brought about by the process of compactive dewatering, with no significant changes due to other processes (e.g. chemical mineralogic adjustment). While it is probable that the model can be extended to other lithologies, it is suspected that chemical adjustments, responding to the effects of time, temperature and pressure on the bulk chemistry of each lithology, will significantly affect changes to xcex and xcexc. Behaviour of the model will also break down if clay dewatering or other processes result in development of an abnormal pore pressure regime in the shale.
Real-world lithologic sequences generally comprise suites of genetically related lithotypes, deposited in a relatively regular manner. Systematic lithologic interpretation from wireline data, using a neural network approach, has been demonstrated (Westphal, Hildegard, and Bornholt, Stephan, (1996, Lithofacies prediction from wireline logs with genetic algorithms and neural networks: Zeitschrift der Deutschen Geologischen Gesellschaft, v. 147, no. 4, p. 465-474), Westphal, Hildegard and Aigner, Thomas, (1997, Seismic stratigraphy and subsidence analysis in the Barrow-Dampier subbasin, northwest Australia: AAPG Bulletin, v.81, No. 10. P. 1721-1749). If volumetric proportions and compactive behavior of each lithotype are known, the contribution of each lithotype to the formation acoustic performance might be computed.             i      .      e      .              xe2x80x83            ⁢      Δ        ⁢          xe2x80x83        ⁢          T      Fmn        =            ∑              i        =        1            n        ⁢                  V        i            *      Δ      ⁢              xe2x80x83            ⁢              T        i            
where xcex94TFmn is the formation ITT (mixed lithologies), Vi is the volumetric proportion of lithotype i, and xcex94Ti is the p wave ITT for lithotype i.
Within reservoir lithologies, porosity variations, which may or may not be systematic, are routine. Detailed analysis of lithotype mixes, and computation of p-wave signal contributions from each lithology, can be expected to only yield a relatively crude estimate of gross interval performance, particularly in complex reservoir lithologies, since only gross lithologic variations, not porosity variations, are defined by the present neural net approach.
Gassman (1951) established p-wave velocity dependence on porosity and density in clastic reservoir lithologies; however, a number of the inputs required to solve the Gassman equation are not routinely available with the precision required for our purposes (rock skeleton bulk and shear modulii, rock grain bulk modulus, interstitial fluid bulk modulus), so an alternative approach to adjustment of computed xcex94T for porosity variations is used.
In most cases, an independent though indirect measure of porosity is available from the resistivity log. If pore fluid properties do not change markedly over a reasonable interval (as is generally the case within genetically related sediment packages not containing hydrocarbons), then the resistivity is a sensitive indicator of porosity variations and is described by the Archie relationship. (Archie, G. E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics, Journal Petroleum Technology 5(1) 54-62). If we use a resistivity device measuring over approximately the same interval as the sonic log, the resistivity data can be transformed to provide a short-interval porosity based modifier to the above composite compactive xcex94TFmn:       Δ    ⁢          xe2x80x83        ⁢          T      Fmn        =                    ∑                  i          =          1                n            ⁢              (                              V            i                    *                      T            i                          )              +          F      ⁢              (        R        )            
where Vi is lithotype volume xcex94Ti is theoretical xcex94T from equation (4) and F(R) is a function (transform) of resistivity.
The presumption here is that for mixed lithologies (mixed sufficiently finely that over the interval of investigation by the sonic tool, discrete layers are not resolved), each lithotype will have proportionate effect on the efficiency of p-wave signal propagation within the interval.
Because we generally lack sufficiently detailed data on the proportions of lithologic components present, this approach is impractical.
Instead, we simplify and generalise the model by regarding any mixed lithology as consisting of shale within which other lithotypes may be mixed or interbedded. Since both rock acoustic performance and electrical resistivity can be generally related to rock porosity (which is itself a function of rock density) we can write:
xcex94TFmn=F(xcfx81) 
and
Rt=Fxe2x80x2(xcfx81) 
Hence in general
xcex94TFmn=Fxe2x80x3(Rt) 
Within pure shale
xcex94TShale=Fxe2x80x3(RShale) 
Taking differences and rearranging
xcex94TFmn=xcex94TShale+[Fxe2x80x3(Rt)xe2x88x92Fxe2x80x3(RShale)]xe2x80x83xe2x80x83(5) 
When applying this equation, we compute xcex94TShale using the compaction based method detailed in equation (4) above, and compute the resistivity difference modifier [Fxe2x80x3(Rt)xe2x88x92Fxe2x80x3(RShale)] using well established petrophysical techniques. This method will now be described.
Formation Resistivity (Rt)
It is presumed that formation resistivity is available as a continuous depth function (a resistivity log).
Shale Resistivity (RShale)
Using general geologic interpretive principles, most probable intervals of cleanest (most nearly xe2x80x9cpurexe2x80x9d) shale are identified on the resistivity log.
Shale resistivity is considered equivalent to the formation resistivity within these intervals. A continuous profile of shale resistivity is developed by interpolation between these intervals. In this manner an expected shale resistivity is established as a continuous function, whether or not any pure shale is actually present.
Functional relationships
The form of resistivity term Fxe2x80x3(Rt) is best established by analysis of Sonic-Porosity and Resistivity-Porosity relationships, and by analysis of observational data compared to computed data.
By analysis of relationships:
From Archie,                               φ          m                =                              aS                          -              n                                ⁢                                    R              w                                      R              t                                                          (        6        )            
where xcfx86 is porosity, m, n and a are constants, S is saturation, Rw and Rt are respectively the resistivity of formation water and of formation.
The general presumption here is that the sequence is fully water saturated i.e. S=1, (see above), so the Archie equation reduces to:       φ    m    =      a    ⁢                  R        w                    R        t            
In equation (6) constants m and a are cementation and tortuosity factors, which have been the subject of considerable previous study, so their ranges and probable values may be estimated with reasonable confidence. Rt is formation resistivity and Rw is formation water resistivity, established by normal log analysis procedures.
But porosity may be directly computed from the sonic log, using either the Wyllie relationship                               φ          sonic                =                              (                                                            Δ                  ⁢                                      xe2x80x83                                    ⁢                  T                                -                                  Δ                  ⁢                                      xe2x80x83                                    ⁢                                      T                    ma                                                                                                Δ                  ⁢                                      xe2x80x83                                    ⁢                                      T                    fl                                                  -                                  Δ                  ⁢                                      xe2x80x83                                    ⁢                                      T                    ma                                                                        )                    ·                      1            C                                              (        7        )            
(wherein xcfx86sonic is porosity computed from sonic ITT. C is a constantxe2x80x94typical range from 1.0 to 1.3 (Raymer, Hunt and Gardner, 1980))
or the Raymer-Hunt-Gardner transform:                               φ          sonic                =                              C            xe2x80x2                    *                      (                                                            Δ                  ⁢                                      xe2x80x83                                    ⁢                  T                                -                                  Δ                  ⁢                                      xe2x80x83                                    ⁢                                      T                    ma                                                                              Δ                ⁢                                  xe2x80x83                                ⁢                T                                      )                                              (        8        )            
where Cxe2x80x2 is a constant (typical value 0.625 to 0.7 according to Alberty (1996) (The influence of the borehole environment upon compression sonic logs; The Log Analyst 37(4) 30-44)).
or the Raiga-Clemenceau equation:                               φ          Sonic                =                  1          -                                    [                                                                    Δ                    ⁢                    T                                    ma                                                  Δ                  ⁢                  T                                            ]                                      1              x                                                          (        9        )            
x typically ranges between 2.3 and 2.4 (Issler, 1992).
We now equate the porosity from the Archie equation to that derived from the various sonic transforms (equations 7, 8 and 9), and solve to express xcex94T in terms of Rt.
Wyllie model:             [                        aS                      -            n                          ⁢                              R            w                                R            t                              ]              1      m        =            [                                    Δ            ⁢                          xe2x80x83                        ⁢                          T              Fmn                                -                      Δ            ⁢                          xe2x80x83                        ⁢                          T              ma                                                            Δ            ⁢                          xe2x80x83                        ⁢                          T              fl                                -                      Δ            ⁢                          xe2x80x83                        ⁢                          T              ma                                          ]        ·          1      C      
from which follows:                     F        xe2x80x3            ⁢              (                  R          t                )              =                  Δ        ⁢                  xe2x80x83                ⁢                  T          Fmn                    =                        Δ          ⁢                      xe2x80x83                    ⁢                      T            ma                          +                              C            ·                          (                                                Δ                  ⁢                                      xe2x80x83                                    ⁢                                      T                    fl                                                  -                                  Δ                  ⁢                                      xe2x80x83                                    ⁢                                      T                    ma                                                              )                        ·                                          (                                                      aS                                          -                      n                                                        ⁢                                                            R                      w                                                              R                      t                                                                      )                                            1                m                                              ⁢                      xe2x80x83                    ⁢          and                                        F        xe2x80x3            ⁡              (                  R          Shale                )              =                  Δ        ⁢                  xe2x80x83                ⁢                  T          Fmn                    =                        Δ          ⁢                      xe2x80x83                    ⁢                      T            ma                          +                  C          ·                      (                                          Δ                ⁢                                  xe2x80x83                                ⁢                                  T                  fl                                            -                              Δ                ⁢                                  xe2x80x83                                ⁢                                  T                  ma                                                      )                    ·                                    (                                                aS                                      -                    n                                                  ⁢                                                      R                    w                                                        R                    shale                                                              )                                      1              m                                          
We substitute these two terms into the above mixed expression for formation xcex94T (equation (5)) to reach the following:                               Δ          ⁢                      xe2x80x83                    ⁢                      T            Fmn                          =                              Δ            ⁢                          xe2x80x83                        ⁢                          T              Shale                                +                      C            ·                          (                                                Δ                  ⁢                                      xe2x80x83                                    ⁢                                      T                    fl                                                  -                                  Δ                  ⁢                                      xe2x80x83                                    ⁢                                      T                    ma                                                              )                        ·                                          (                                                      aS                                          -                      n                                                        ⁢                                      R                    w                                                  )                                            1                m                                      ·                          [                                                                    (                                          1                                              R                        t                                                              )                                                        1                    m                                                  -                                                      (                                          1                                              R                        Shale                                                              )                                                        1                    m                                                              ]                                                          (        10        )            
Using similar logic, we develop from the Raymer-Hunt Gardner model                               Δ          ⁢                      xe2x80x83                    ⁢                      T            Fmn                          =                              Δ            ⁢                          xe2x80x83                        ⁢                          T              Shale                                +                      Δ            ⁢                          xe2x80x83                        ⁢                                          T                ma                            ·                              [                                                      1                                          (                                              1                        -                                                                              1                                                          C                              xe2x80x2                                                                                ⁢                                                                                    (                                                                                                aS                                                                      -                                    n                                                                                                  ⁢                                                                                                      R                                    w                                                                                                        R                                    t                                                                                                                              )                                                                                      1                              m                                                                                                                          )                                                        -                                      1                                          (                                              1                        -                                                                              1                                                          C                              xe2x80x2                                                                                ⁢                                                                                    (                                                                                                aS                                                                      -                                    n                                                                                                  ⁢                                                                                                      R                                    w                                                                                                        R                                    Shale                                                                                                                              )                                                                                      1                              m                                                                                                                          )                                                                      ]                                                                        (        11        )            
and from the Raiga-Clemenceau model:                               Δ          ⁢                      xe2x80x83                    ⁢                      T            Fmn                          =                              Δ            ⁢                          xe2x80x83                        ⁢                          T              Shale                                +                      Δ            ⁢                          xe2x80x83                        ⁢                                          T                ma                            ·                              [                                                      1                                                                  (                                                  1                          -                                                                                    (                                                                                                aS                                                                      -                                    n                                                                                                  ⁢                                                                                                      R                                    w                                                                                                        R                                    t                                                                                                                              )                                                                                      1                              m                                                                                                      )                                            x                                                        -                                      1                                                                  (                                                  1                          -                                                                                    (                                                                                                aS                                                                      -                                    n                                                                                                  ⁢                                                                                                      R                                    w                                                                                                        R                                    Shale                                                                                                                              )                                                                                      1                              m                                                                                                      )                                            x                                                                      ]                                                                        (        12        )            
Hence the inventor finds that the general form of a synthetic sonic algorithm is:
xcex94TFmn=xcex94TShale+[Fxe2x80x3(Rt)xe2x88x92Fxe2x80x3(RShale)]xe2x80x83xe2x80x83(5) 
where the resistivity terms are computed in the manner of equations (10), (11), and (12) above, and xcex94TShale is computed from the compactive sonic theoretical algorithm defined above in equation (4).
In one form, although it need not be the only or indeed the broadest form, the invention resides in a method of computing a theoretical sonic log including the steps of:
computing an ideal theoretical shale sonic log xcex94TShale; and
correcting the ideal theoretical shale sonic log with measured resistivity data using the relation xcex94TFmn=xcex94TShale+F(Res) where F(Res) is a selected resistivity function.
In a further form the method includes the steps of:
calculating a compactive constant k;
estimating a modulus function (xcex+2xcexc) as a function of depth using the constant k and reliable observations of interval transit time (xcex94T);
calculating the shale interval transit time xcex94TShale;
calculating a resistivity modifier term which is then added to the shale interval transit time
The compactive constant k is suitably determined by using point density data to solve
xcfx81(z)=xcfx81∞+(xcfx810xe2x88x92xcfx81∞)exe2x88x92kz 
The point density data xcfx81(z) may be derived from core measurements or log measurements.
The step of estimating the modulus function (xcex+2xcexc) is suitably performed by solving
xcex94T(z)2*(xcexz+2xcexcz)=(A*xcex94T∞2)+(B*xcex94T02xe2x88x92A*xcex94T∞2)exe2x88x92kz 
by substituting with observed xcex94T(z) at various z and interpolating between calculations.
The step of calculating the shale interval transit time xcex94Tsh is performed by solving       Δ    ⁢          xe2x80x83        ⁢          T      ⁡              (        z        )              =            [                        [                                    (                              A                *                Δ                ⁢                                  xe2x80x83                                ⁢                                  T                  ∞                  2                                            )                        +                                          (                                                      B                    *                    Δ                    ⁢                                          xe2x80x83                                        ⁢                                          T                      o                      2                                                        -                                      A                    *                    Δ                    ⁢                                          xe2x80x83                                        ⁢                                          T                      ∞                      2                                                                      )                            ⁢                              ⅇ                Ks                                              ]                          (                                    λ              s                        +                          2              ⁢                              μ                s                                              )                    ]              1      2      
The step of determining shale resistivity is suitably performed by taking values of shale resistivity from a resistivity log, and interpolating between points of observation to yield a continuous estimate of shale resistivity.
Suitably the interval transit time is corrected using the transforms of resistivity defined in equations 10, 11 and 12 above.
In another form, the invention resides in a method of calibrating an acquired sonic log including the steps of:
computing an ideal theoretical sonic log xcex94TShale;
comparing the ideal theoretical sonic log to high confidence regions of the acquired sonic log;
modifying the ideal theoretical sonic log to correlate with said high confidence regions; and
substituting the ideal theoretical sonic log for the acquired sonic log in regions other than the high confidence regions.